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Abstract. CG, SYMMLQ, and MINRES are Krylov subspace methods for solving symmetric 
systems of linear equations. When these methods are applied to an incompatible system (that is, 
a singular symmetric least-squares problem), CG could break down and SYMMLQ's solution could 
explode, while MINRES would give a least-squares solution but not necessarily the minimum-length 
(pseudoinverse) solution. This understanding motivates us to design a MINRES-like algorithm to 
compute minimum-length solutions to singular symmetric systems. 

MINRES uses QR factors of the tridiagonal matrix from the Lanczos process (where R is upper- 
tridiagonal). MINRES-QLP uses a QLP decomposition (where rotations on the right reduce R to 
lowcr-tridiagonal form). On ill-conditioned systems (singular or not), MINRES-QLP can give more 
accurate solutions than MINRES. We derive preconditioned MINRES-QLP, new stopping rules, and 
better estimates of the solution and residual norms, the matrix norm, and the condition number. 
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1. Introduction. We are concerned with iterative methods for solving a sym- 
metric linear system Ax = b or the related least-squares (LS) problem 

min||a;||2 s.t. x € argmin || A.t — 6||2, (1-1) 

X 

where A E M"^" is symmetric and possibly singular, b G M", A ^ 0, and b ^ Q. Most 
of the results in our discussion are directly extendable to problems with complex 
Hermitian matrices A and complex vectors b. 

The solution of (1-1), called the minimum-length or pseudoinverse solution [18], 
is formally given by = {A'^A)'' A^b = {A'^)'< Ab = {A'i)^Ab, where denotes the 
pseudoinverse of A. The pseudoinverse is continuous under perturbations E for which 
rank {A + E) = rank (A) [49], and a;^ is continuous under the same condition. Problem 
(1.1) is then well-posed [19]. 

Let A = UMJ^ be an eigenvalue decomposition of A, with U orthogonal and 
A = diag(Ai, . . . , A„). We define the condition number of A to be k{A) = nii™^^^ojA | ' 
and we say that A is ill-conditioned if n{A) 3> 1. Hence a singular matrix co'uld he 
well-conditioned or ill-conditioned. 
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SYMMLQ and MINRES [39] are Krylov subspace methods for solving symmetric 
indefinite systems Ax = b. SYMMLQ is reliable on compatible systems even if A is 
ill-conditioned or singular, while on (singular) incompatible problems its iterates Xk 
diverge to a multiple of a nuUvector of A [10, Proposition 2.15] and [10, Lemma 2.17]. 
MINRES seems more desirable to users because its residual norms are monotonically 
decreasing. On singular compatible systems, MINRES returns x^ (see Theorem 3.1). 
On singular incompatible systems, MINRES is reliable if terminated with a suitable 
stopping rule involving ||v4rfe|| (see Lemma 3.3), but the solution will not be x^ . 

Here we develop a new solver of this type named MINRES-QLP [10]. The aim 
is to deal reliably with compatible or incompatible systems and to return the unique 
solution of (1.1). We give theoretical reasons why MINRES-QLP improves the accuracy 
of MINRES on ill-conditioned systems, and illustrate with numerical examples. 

Incompatible symmetric systems could arise from discretized semidefinite Neu- 
mann boundary value problems [27, section 4], and from any other singular systems 
involving measurement errors in b. Another potential application is large symmetric 
indefinite low-rank Toeplitz LS problems as described in [16, section 4.1]. 

1.1. Notation. The letters z, j, k denote integer indices, c and s cosine and sine 
of some angle 9, the fcth unit vector, e a vector of all ones, and other lower-case 
letters such as 6, m, and x (possibly with integer subscripts) denote column vectors. 
Upper-case letters A, T^, Vk, ■ ■ ■ denote matrices, and Ik is the identity matrix of 
order k. Lower-case Greek letters denote scalars; in particular, e sa 10"^^ denotes the 
floating-point precision. If a quantity 6k is modified one or more times, we denote its 
values hy 6k, Sf, , Sj^ , ... . The symbol || • || denotes the 2-norm of a vector or matrix. 
For an incompatible system. Ax w 6 is shorthand for the LS problem min^ ||j4x — &|j. 

1.2. Overview. In sections 2-4 we briefly review the Lanczos process, MINRES, 
and QLP decomposition before introducing MINRES-QLP in section 5. We derive 
norm estimates in section 6 and preconditioned MINRES-QLP in section 7. Numerical 
experiments are described in section 8. 

2. The Lanczos process. Given A and 6, the Lanczos process [30] computes 
vectors Vk and tridiagonal matrices Tk according to Vq = 0, = b, and then^ 

Pk = Avk, ak = WfePfc, Pk+iVk+i ^Pk- ctkVk - PkVk^i 
for fc = 1, 2, . . . where we choose > to give \\vk\\ = 1. In matrix form. 



AVk = Vk+iTk, Tk = 



In exact arithmetic, the columns of 14 are orthonormal and the process stops with 
k — i and Pe+i = for some i < n, and then AVi = VgTg. For derivation purposes 
we assume that this happens, though in practice it is unlikely unless 14 is reorthog- 
onalized for each k. In any case, (2.1) holds to machine precision and the computed 
vectors satisfy l|Vfc]|i sa 1 (even ii k ^ n). 



"1 



/92 
2 Ol2 



&4-1J 



Tk 



, Vk=[v^ 



Vk\ 



(2.1) 



1 Numerically, = Av^ - Pk^k-l, "ft = ''fePft. Pk+l'Uk+l = Pk - ^kVk is slightly better [38]. 
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2.1. Properties of the Lanczos process. The fcth Krylov subspace generated 
by A and b is defined to be ICk{A, b) = span{6, Ab, A^b, A'^^^b} = span(Vfe). The 
following properties should be kept in mind: 

1. If ^ is changed to A — al for some scalar shift a, Tk becomes Tk—al and Vk is 
unaltered, showing that singular systems are commonplace. Shifted problems 
appear in inverse iteration or Rayleigh quotient iteration. 

2. Tk has full column rank k for all k < £. 

3. If A is indefinite, some Tk might be singular for k < £, but then by the Sturm 
sequence property (see [18]), Tk has exactly one zero eigenvalue and the strict 
interlacing property implies that Tk±i are nonsingular. Hence Tk cannot be 
singular twice in a row (whether A is singular or not). 

4. Ti is nonsingular if and only if & € range(A). (See appendix A.) 

3. MINRES. Algorithm MINRES [39] is a natural way of using the Lanczos 
process to solve Ax = b or min^ \\Ax — b\\ . For k < i, if Xk = VkUk for some vector yk, 
the associated residual is 

Tk^b - Axk = & - AVkUk = - Vk+iTkVk = Vk+iiPiCi - Tkjjk)- (3.1) 

To make small, it is clear that j3 iCi ~ Tkyk should be small. At this iteration k, 
MINRES minimizes the residual subject to Xk € /Cfc(A, b) by choosing 

= arg min ||Tfc?/ - /3iei||. (3.2) 



This subproblem is processed by the expanding QR factorization: Qq = 1 and 



ik-i 

Ck Sfc 
Sfc — Cfc 



{k,k+l 



1 



QklTk /3iei 



^fe tk 

0fe 



, (3.3) 



where Ck and Sk form the Householder reflector Qk,k+i that annihilates /3fc+i in Tk to 
give upper-tridiagonal Rk, with Rk and tk being unaltered in later iterations. 

When k < £, the unique solution of (3.2) satisfies RkVk = tk- Instead of solving 
for yk , MINRES solves RJId'^ = by forward substitution, obtaining the last column 
dk of Dk at iteration k. At the same time, it updates Xk via a;o = and 

Xk = Vfej/fe = DkRkVk = Dktk = Xk-i + Tkdk, Tk = eltk. (3.4) 

When k = i, we can form Ti but nothing else expands. In place of (3.1) and 
(3.3) we have re = Ve{Piei — Tgye) and Qi-i[Tg Piei] = [Re te] and it is natural 
to choose ye from the subproblem 

min IjT^y^ - /3iei|| = min - . (3.5) 

There are two cases to consider: 

1. If Ti is nonsingular, Reye = te has a unique solution. Since AVeye = VeTeye = 
6, the problem is solved by xe = Veye with residual re = (the system is 
compatible, even if A is singular). Theorem 3.1 proves that xe = xK 

2. If Te is singular, A and Re are singular {Ree = 0) and both Ax — b and 
RiUi — te are incompatible. This case was not handled by MINRES in [39]. 
Theorem 3.2 proves that the MINRES point xe-i is a least-squares solution 
(but not necessarily x^). Theorem 5.1 proves that the MINRES-QLP point 
Xe — Veyl = x\ where y^ is the min-length solution of (3.5). 
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3.1. Further details of MINRES. To describe MINRES-QLP tliorouglily, we 
need furtlier details of the MINRES QR factorization (3.3). For 1 < k < £, 



Rk 





71 



(2) 
72 



£3 



7f' 





n 




Cl 








S1C2 




'tk 






= /3i 


















Tk 




Si • • • Sk-lCk 




_4'k_ 




_Sl • • • Sfc_lS/c_ 



(3.6) 





(where the superscripts are defined in section 1.1). With (j>o = /3i > 0, the full action 
of Qk.k+i in (3.3), including its effect on later columns of Tj, k < j < £, is described by 



Cfe 
Sk 



Sk 



Ik 

Pk+1 



5k+l 
Oik+1 





Pk+2 













(2) 

Tk 



^(2) 

7fc+i 



efc+2 

4+2 



Tk 
4>k 



(3.7) 



(2) 

where Sk ~ (3k+i/\\['yk Pk+i]\\ > 0, giving 71, 7^ > with Rj nonsingular for each 
j < k < i. Thus the dj in (3.4) can be found from 

= t;i/7i, d2 



{v2-62d,)/jp, 

(vj - sfUj^i - ejdj^2)hf\ 



Also, Tfe = 4>k-iCk and 0^ = 4>k-iSk > 0. Hence from (3.1)-(3.3), 

Ik/cll = ||Tfc2/fe - /3iei|| = (/)fe ||rfe|| = ||rfe_i||sfc, 



(3.8) 



(3.9) 



which is nonincreasing and tending to zero if Ax = b is compatible. 

Remark 3.1. If k < £ and is singular, we have 7^ = 0, Sk — 1, and ||rfc|| = 
||rfc_i|| (not a strict decrease), but this cannot happen twice in a row (cf. section 2.1). 

Remark 3.2. IfTgis singular, MINRES sets the last element of ye to be zero. The 
final point and residual stay as xi-i andre^i with\\ri^i\\ — (f>i-i — j3isi ■ ■ ■ si-i > 0. 

3.2. Compatible systems. The following theorem assures us that MINRES is 
a useful solver for compatible linear systems even if A is singular. 

Theorem 3.1 ([10, Theorem 2.25]). Ifb€ range(A), the final MINRES point xg 
is the minimum-length solution of Ax = b (and rg = b — Axi — 0). 

Proof. If 6 e range(yl), the Lanczos process gives AVe — VfTi with nonsingular 
Ti, and MINRES terminates with Axi = b and xi = Viyi = Aq, where q = VfT^^yi. 
If some other point x satisfies Ax = 6, let p = x ~ xi. We have Ap = and 



T 

xjp 



q^Ap = Q. Hence Ijirlp = + p|| 



\xe\ 



X ~ Xi. 

f 2xjp H 



\\pr>\\xer. □ 



3.3. Incompatible systems. For a singular LS problem Ax ~ b, the optimal 
residual vector r is unique, but infinitely many solutions x give that residual. In the 
following example, MINRES finds a least-squares solution (with optimal residual) but 
not the minimum-length solution. 

Example 3.1. Let A = diag(l, 1,0) and b = e. The minimum-length solution to 
Ax K, b is x^ — [1 1 0]^ with residual r = b — Ax^ = 63 and Ar — 0. MINRES returns 
the solution x^ — e (with residual r^ ^ b — Ax"^ = e^ ~ r and Ar^ = 0)- 

Theorem 3.2 ([10, Theorem 2.27]). Ifb^ range(A), then \\Are^i\\ = and the 
MINRES xi_i is an LS solution (but not necessarily x^ ). 

Proof. Since b ^ range(A), Ti is singular and Ru = — 0. By Lemma 3.3 below, 
A(Ax{^i — 6) = —Are^i = —\\re-i\\"feV£ = 0. Thus xg^i is an LS solution. □ 
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3.4. Norm estimates in MINRES. For incompatible systems, rk (3.1) will 
never be zero. However, all LS solutions satisfy A^x — Ab, so that Ar = 0. We 
therefore need a new stopping condition based on the size of ||Arfe||. In applications 
requiring nullvectors, || Aa;fc|| is also useful. We present efficient recurrence relations for 
||Arfc|| and ||^a;fe|| in the following Lemma, which was not considered in the framework 
of MINRES when it was originally designed for nonsingular systems [39] . 

Lemma 3.3 (Ark, \\Ark\\, \\Axk\\ for MINRES). Ifk<i, 

Ark Ikfell ilk+iVk+i + 4+2W/C+2) {where Sk+2Vk+2 = ijk = l-l), 

i,l = \\Ark\\^ - \\rkf {hk+i? + [5k+2?) (where 4+2 = i/fc = ^ - 1), 
ivl = \\Axk\\^ ^ujl_i+Ti, uJo = 0. 

Proof. For k < Rk is nonsingular. From (3.1)-(3.4) with Rkyk = tk we have 

Vk^ = (pkVk+iQlck+i, (3.10) 

Ark = 4>kVk+2Tk±iQlek+i, 

Qk Tk+l ^ — Qk [Tk+l l3k+2ek+l] = Qk 

el+iQkTk+i^ ^ [0 7fc+i 4+2], 
see (3.7), where AVk+i = Vk+iTk+i and we take 6k+2 = if fc = £ — 1, so 

Ark ^ (f>kVk+2[0 7fc+i Sk+2] ^ (pk (jk+iVk+i + Sk+2Vk+2) , 
i;l^\\Ark\\' = \\rkr{hk+ir + [Sk+2]'). 

For the recurrence relations of Axk and its norm, we have 

Axk = AVkVk = Vk+iTkiik = Vk+iQl 
ujI = \\Axkf = \\tkf = \\tk-if +tI=ujI_^+tI 

□ 

Note that even using finite precision the expression for ■0^ is extremely accurate 
for the versions of the Lanczos algorithm given in section 2, since (taking \\vj\\ — 1 
with negligible error), WArkW^ = (j)l([^k+if + '2.lk+i5k+2vl^^Vk+2 + [4+2]^), where 
from (3.7) \5k+2\ < Pk+2, while from [38, (18)] l3k+2\vlj^iVk+2\ < 0(£)\\A\\, and with 
|7fc+i| < ||A||, see [38, (19)], we see that \ik+i5k+2vl+^Vk+2\ < 0(e)||A||2. 

Typically H^J'fcll is not monotonic, while clearly ||rfe|| and ||Axfe|| are monotonic. 
In the eigensystem A — UAU'^, let U — [Ui U2], where the eigenvectors Ui cor- 
respond to nonzero eigenvalues. Then Pa = UiUi and P^ = U2U2 are orthogonal 
projectors [53] onto the range and nuUspace of A. For general linear LS problems, 
Chang et al. [7] characterize the dynamics of \\rk\\ and in three phases defined 

in terms of the ratios among ||rfc||, ||PA^fc||j and ll^'j^^J'fcll, and propose two new stop- 
ping criteria for iterative solvers. The expositions [1, 26] show that these estimates 
are cheaply computable in CGLS and LSQR [40, 41]. 



Tk = Vk+iQl' 
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3.5. Effects of rounding errors in MINRES. MINRES should stop if Rk is 
singular (which theoretically implies k = £ and A is singular). Singularity was not 
discussed by Paige and Saunders [39], but they did raise the question: Is MINRES 
stable when Rk is ill-conditioned? Their concern was that \\Dk\\ could be large in 
(3.8), and there could then be cancellation in forming Xk-i + T^dk in (3.4). 

Sleijpen, Van der Vorst, and Modersitzki [47] analyzed the effects of rounding 
errors in MINRES and reported examples of apparent failure with a matrix of the 
form A = QDQ^, where D is an ill-conditioned diagonal matrix and Q involves a 
single plane rotation. We were unable to reproduce MINRES's performance on the 
two examples defined in Figure 4 of their paper, but we modified the examples by using 
an n X n Householder transformation for Q, and then observed similar difficulties with 
MINRES — see Figure 8.2. The recurred residual norm (f)^^ is a good approximation 



until the last few iterations. 



The recurred norms (pjf 



to the directly computed 

then keep decreasing but the directly computed norms ||r^^' || become stagnant or even 
increase (see the lower subplots in Figure 8.2). 

Remark 3.3. Note that we do want (pk to keep decreasing on compatible systems, 
so that the test (j)^ < toZ(|| A|| ||a;j.|| + ||&||) with tol > e will eventually be satisfied even 
if the computed \\rk\\ is no longer as small as 4>k- 

The analysis in [47] focuses on the rounding errors involved in the n lower tri- 
angular solves RJ^DJ^ ~ (one solve for each row of Dk), compared to the single 
upper triangular solve RkVk = tk (followed by Xk = Vkyu) that would be possible at 
the final k if all of Vk were stored as in GMRES [44] . We shall see that a key feature 
of MINRES-QLP is that a single lower triangular solve suffices with no need to store 
Vfc, much the same as in SYMMLQ. 



4. Orthogonal decompositions for singular matrices. In 1999 Stewart pro- 
posed the pivoted QLP decomposition [51], which is equivalent to two consecutive QR 
factorizations with column interchanges, first on A, then on R^: 



QrAUr - 



'R 


s 




-rt 


0' 


Hi = 


'R 


0" 








Ql 














(4.1) 



giving nonnegative diagonal elements, where 11/? and 11^ are permutations chosen to 
maximize the next diagonal element of R and R at each stage. This gives A — QLP, 
where 



L = 



R"^ 




with Q and P orthogonal. Stewart demonstrated that the diagonals of L (the L- 
values) give better singular- value estimates than the diagonals of R (the R-values), 
and the accuracy is particularly good for the extreme singular values ai and cr„: 



Ri. 



CTj, fi > maxLii > maxRii, minRa > imnLa > an- (4.2) 



The first permutation II f( in pivoted QLP is important. The main purpose of the 
second permutation 11]^ is to ensure that the L- values present themselves in decreasing 
order, which is not always necessary. If Hji = 11^ = /, it is simply called the QLP 
decomposition. 
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5. MINRES-QLP. We now develop MINRES-QLP for solving ill-conditioned or 
singular symmetric systems Ax w b. The Lanczos framework is the same as in MIN- 
RES, but we handle Ti in (3.5) with extra care when it is rank-deficient. In this case, 
the normal approach to solving (3.5) is via a QLP decomposition of Ti to obtain the 
(unique) minimum-length solution yi [51, 18]. Thus in MINRES-QLP we use a QLP 
decomposition of Tk in subproblem (3.2) for all k < £. This is the MINRES QR (3.3) 
followed by an LQ factorization of Rk'- 



QkTk — 



Rk 




RkPk = Lk 



so that QkTkPk = 



Lk 




(5.1) 



where Qfc and are orthogonal, R^. is upper tridiagonal and L^. is lower tridiagonal. 
When k < £, Rk and Lk are nonsingular. MINRES-QLP obtains the same solution as 
MINRES, but by a different process (and with different rounding errors). Defining u 
by y = Pku, we see from (3.3) that 



Qk{Tky - Piei) 



Lk 




and (3.2) is solved by L^Uk ~ and j/j. = PkU^- The MINRES-QLP estimate of x is 
therefore Xk = V^yk = VkPk^k — W^Uk, with theoretically orthonormal Wk = V^Pk- 
We will see that only the last three columns of Vk are needed to update Xk- 

5.1. The QLP factorization of Tk. The QLP decomposition of each Tk must 
be without permutations in order to ensure inexpensive updating of the factors as k 
increases. Our experience is that the desired rank-revealing properties (4.2) tend to 
be retained, perhaps because of the tridiagonal structure of Tk and the convergence 
properties of the underlying Lanczos process. 

For k < £, the QLP decomposition of Tk (5.1) gives nonsingular tridiagonal Rk 
and Lk- As in MINRES, Qk is a product of Householder reflectors, see (3.3) and (3.7), 
while Pk involves a product of pairs of essentially 2x2 reflectors: 



{k,k + l 



Q3,4 Q2,3 Ql,2, 



Pk = P^ 



1,2 



Pi.:iP. 



2,3 



Pk~2,kPk-l.k- 



For MINRES-QLP to be efficient, in the fcth iteration {k > 3) the application of the 
left reflector Qk.k+i is followed immediately by the right reflectors Pk-2.k, Pk-i,k, so 
that only the last 2x2 principal submatrix of the transformed Tk will be changed 
in future iterations. These ideas can be understood more easily from Figure 5.1 and 
the following compact form, which represents the actions of right reflectors on Tk 
(additional to Qk,k+i (3.7)): 



■ (5) 

^k-l 
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Sk3 
~Ck3. 



(5) 

Tk-i 



^k 



lk 



(4) 



(5.2) 



5.2. The diagonals of Lk- Figure 5.2 shows the relation between the singular 
values of A and the diagonal elements of Rk and Lk with k = 19. This illustrates 
(4.2) for matrix ID 1177 from [54] with n = 25. 
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Fig. 5.1. QLP with left and right reflectors interleaved on T5. This figure can be reproduced 
with the help of QLPfigS.m. 



nonzero a{A) = \X{A)\, n = 25 7^\fc = 19 




7f 



20 





Fig. 5.2. Upper left: Nonzero singular values of A sorted in decreasing order. Upper middle 
and right: The diagonals 7^^ of (red circles) from MINRES are plotted as red circles above or be- 
low the nearest singular value of A. They approximate the extreme nonzero singular values of A well. 
Lower: The diagonals 7^ of L^^ (red circles) from MINRES-QLP approximate the extreme nonzero 
singular values of A even better. An implication is that the ratio of the largest and smallest diagonals 
of Lk provides a good estimate ofK{A). To reproduce this figure, run test_minresqlp_fig3(2). 
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5.3. Solving the subproblem. With yk ~ PkUk, subproblem (3.2) becomes 

(5.3) 



Uk 



arg mm 





Lk 


u — 


tk 
















where tk and (pk are as in (3.3) and (3.6). At the start of iteration k, the first k — 3 
elements of Uk, denoted by fij for j < k~3, are known from previous iterations; see 
the 10th matrix in Figure 5.1. The remainder depend on the rank of Lk- 

1. If rank(Lfe) — k (so k < £, or k — £ and b G range(A)), we need to solve the 
last three equations of LkUk — tk'. 
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Tk-1 



(5.4) 

If fc = ^ and b ^ range(A), the last row and column of Lk are zero, and we 
only need to solve the last two equations of Lk-iUk-i = tk-i, where 



Lk 



Lk-i 




Uk-l 







, Uk = 





1 



(6) 

rk-2 

^(2) ^(5) 





" (3) ■ 








A^fc-2 




Tk-2 




(2) 




fk-l_ 



The corresponding solution estimate is Xk = VkVk — ^kPkUk = WkUk, where 
Wk = VkPk = [Vk-iPk-i Vk\Pk-2,kPk-i,k 



(4) 



fc-3 



(3) 
%-2 



W 



(2) 
k-1 



Vk 



Pk-2,kPk-l,k 



(4) (3) (2) 

^k-2 ^k-1 ^k 



W^Wk = Ik, rangeiWk) = ICk{A, b), 
and we update Xk-2 and compute Xk by short-recurrence orthogonal steps: 



(2) 
''fe-2 



Xk 



(2) 
''fe-3 

(2) 
''fe-2 



^k-2^^k-2^ where x 



(2) 
fc-3 



^'^fe-3"fe-3' 



,,,(3) (2) 
^k-lf^k-1 



(2) 



(5.5) 
(5.6) 

(5.7) 

(5.8) 
(5.9) 



5.4. Termination. When k — i, Qk.k+i is not formed or applied, see (3.3) 
and (3.7), and the QR factorization stops. In MINRES-QLP, we still need to apply 
Pk-2,kPk-i,k on the right to obtain the minimum-length solution; see Figure 5.1. 



Theorem 5.1 ([10, Theorem 3.1]). In MINRES-QLP, xi 



Proof. When b G range(^), the proof is the same as that for Theorem 3.1. 

When 6 i range(v4), for aU u = [ui-i Hkf G M*^ that solves (5.3), MINRES-QLP 
returns the min-length LS solution = [ug-i 0]^by the construction in (5.5). For 
any x € Ta,nge{Wi) = JCe{A, b) by (5.7), 



\Ax-b\\ 



AWeu-b\\ 
Qe-iTiPiU 



WAVePiU - b\\ = WVeTePiU - f3iVtei\\ = \\TtPiu - pici 



te-i 



Le-i 




te- 



Since Lg^i is nonsingular, (j)e-i = min||j4a; — b\\ can be achieved by xi = Weug = 
Wg-iui^i and \\xi\\ = || Wf-iUf-i || by (5.7). Thus Xf is the min-length LS 
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solution of \\Ax — b\\ in JCe{A,b), i.e., xg ~ argniin{||a;|| | A^x = Ab, x e JC£{A,b)}. 
Likewise yi = Pgui is the min-length LS solution of /3iei|| and so yi € range(rf), 

i.e. yi = Tiz for some z. Thus Xi — Viyi = VfTiz = AVgz G range(^). We know 
that x^ = argminlllxll | A^x = Ab, x £ M"} is unique and x^ G range(A). Since 
xi € range(A), we must have xe ^ xK □ 

5.5. Transfer from MINRES to MINRES-QLP. On well-conditioned sys- 
tems, MINRES and MINRES-QLP behave very similarly. However, MINRES-QLP re- 
quires one more vector of storage, and each iteration needs 4 more axpy's (y -s— ax + y) 
and 3 more vector scalings {x ^ ax). Thus it would be a desirable feature to invoke 
MINRES-QLP from MINRES only if A is ill-conditioned or singular. The key idea is to 
transfer to MINRES-QLP at an iteration where Tk is not yet too ill-conditioned. The 
MINRES and MINRES-QLP solution estimates are the same, so from (3.4), (5.9), and 
(5.3): xf Xk Dktk = WkUk = WkL^^tk- Now from (3.8), (5.1), and (5.6), 

DkLk = {VkRk^){RkPk) = VkPk = Wk, (5.10) 

and the last three columns of Wk can be obtained from the last three columns of 
Dk and Lk- (Thus, we transfer the three MINRES basis vectors dk-2,dk-i,dk to 
Wk-2,Wk-i,Wk-) In addition, we need to generate Xf._2 using (5.8): 

(2) M (3) (2) (2) 

It is clear from (5.10) that we still need to do the right transformation RkPk = Lk 
in the MINRES phase and keep the last 3x3 principal submatrix of Lk for each k 
so that we are ready to transfer to MINRES-QLP when necessary. We then obtain 
a short recurrence for ||xfe|| (see section 6.5) and for this computation we save flops 
relative to the original MINRES algorithm, where ||xfc|| is computed directly. 

In the implementation, the MINRES iterates transfer to MINRES-QLP iterates 
when an estimate of the condition number of Tk (see (6.3)) exceeds an input param- 
eter trancond. Thus, trancond > 1/e leads to MINRES iterates throughout, while 
trancond = 1 generates MINRES-QLP iterates from the start. 

5.6. Comparison of Lanczos-based solvers. We compare MINRES-QLP with 
CG, SYMMLQ, and MINRES in Tables 5.1-5.2 in terms of subproblem definitions, 
basis, solution estimates, flops and memory. A careful implementation of SYMMLQ 
provides a point in ]Ck+i{A,b) as shown. All solvers need storage for Vk, Vk+i-, Xk, 
and a product pk — Avk each iteration. Some additional work-vectors are needed for 
each method (e.g., dk^i and dk for MINRES, giving 7 work-vectors in total). 

6. Stopping conditions and norm estimates. This section derives several 
norm estimates that are computed in MINRES-QLP. As before, we assume exact 
arithmetic throughout, so that Vk and Qk are orthonormal. Table 6.1 summarizes 
how the norm estimates are used to formulate three groups of stopping conditions. 
The second NRBE test ||v4ri:|| < ||A||||rfe||toZ is from Stewart [50] with symmetric A. 

6.1. Residual and residual norm. First we derive recurrence relations for rk 

and its norm (j)k = WrkW- 

Lemma 6.1 (r^ and ||rfc|| for MINRES-QLP and monotonicity of ||rfc||). 

• Ifk< £, then rank(Lfc) = k, rk ^ slrk-i-cj)kCkVk+i, and (j)k = (t)k-iSk > 0. 

• //rank(_L^) £, then = 0. 

• //rank(Lf) = £ — 1, then rg = r£_i ^ 0, and \\r£\\ = 4>{-i > 0. 
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Table 5.1 

Subprohlems defining for CG, SYMMLQ, MINRES, and MINRES-QLP. 



Method 



Subprobleni 



Factorization Estimate of 



cgLanczos 
[24, 39, 48] 



Cholesky: 



Xk = VkVk 
e/Cfe(A,6) 



SYMMLQ 
[39, 45] 

MINRES 

m 

MINRES-QLP 
[10] 



ijk+i = argminygK'=+i WvW 
s.t. T^y = /3iei 

yk = arg min \\Tk_y ~ (3iei\\ 



LQ: 

Tk^Ql= [Lk 0] 
QR: 



Vk = argmiHygufc ||y|| 

s.t. y e arg min ||7fey - /?iei| 



QkTk 

QLP: 

QkTkPk ~ 



Rk 




Lk 




x^ Vk+iVk+i 
elCk+i{A,b) 

xf = VkVk 
elCk{A,b) 

x1 = VkVk 
elCk{A,b) 



Table 5.2 

Bases, subproblem solutions, storage, and work for each method. 



Method 


New basis 


Zk,tk,Uk 


Xk estimate 


vecs flops 


CgLanczos 


Wk = VkL-^ 


LkDkZk = Piei 


x^ = WkZk 


5 8n 


SYMMLQ 


Wk^Vk+^Ql\^^^ 


LkZk^ Piei 


x^^WkZk 


6 9n 


MINRES 


Dk = VkR^^ 


tk^l3i[lk OjQfcCi 


xf^Dktk 


7 9n 


MINRES-QLP 


Wk = VkPk 


LkUk=l3i[lk OjQfcei 


x^=WkUk 


8 14?! 



Proof, li k < £, the residual is the same as for MINRES. We have ||rfc| 
bk-iSk > 0; see (3.6)-(3.9). Also from = (f>kVk+iQlek+i (3.10) we have 



1 



1 



'4-1 



rk = (f>k \Vk Vk+i 

= 4>k [Vk Vk+i 

slrk-i - 4)kCkVk+i by (3.10) 



Cfc Sfc 

Sfe -Cfc. 





"0" 









_1_ 



SkSk 
-Ck 



\Vk 



Vk+l\ 



by (3.3), 

-Ck 



(t>kSkVkQk-iek - 4>kCkVk+i = 4>k-islVkQl^iek - 4>kCkVk+i 



If Ti is nonsingular, = 0. Otherwise Qf-i/ has made the last row of Ri zero, so the 
last row and column of are zero; see (5.5). Thus ri = r^_i ^ 0; see Remark 3.2. □ 

6.2. Norm of Ark- Next we derive recurrence relations for Ark and its norm 
ifik = ll^rfcll, and we show that Ark is orthogonal to JCk{A,h). 
Lemma 6.2 {Ark and Vfc = ||-4rfc|| for MINRES-QLP). 

• If k < i, then rank(Lfc) = k, Ark = \\rk\\{lk+iVk+i + Sk+2Vk+2) and i/^fe = 
||rfc||||[7fc+i 4+2]||, where 4+2 = if k = £-1. 

• //rank(L^) = £, then Ari = and ipg = 0. 

• //rank(L^) — £^1, then Ar£ — Ar£_i = 0, and WipiW = ipt-i = 0. 

Proof. For the first case, the proof is essentially the same as the proof of Lemma 3.3. 
For the other two cases, the results follow directly from Lemma 6.1. □ 



12 



S.-C. T. CHOI, C. C. PAIGE, AND M. A. SAUNDERS 



Table 6.1 

Stopping conditions in MINRES-QLP. NRBE means normwise relative backward error, and 
tol, maxit, maxcond, and maxxnorm are input parameters. All norms and k{A) are estimated by 
MINRES-QLP. 



Lanczos 


NRBE 


Regularization attempts 


/3fe+i < n\\A\\e 
k = maxit 


||rfe||/(||A||||a;fe|| + ||6||) < tol 
WAvkW/ i\\A\\\\rk\\) < tol 


k(A) > maxcond 
\\xk\\ > maxxnorm 



6.3. Matrix norms. For Lanczos-based algorithms, ||v4|| > HV'^j^j^AVjtH = HTfelj. 
Define 

= 0, A^^'> = max^ |||T,-ej||} = max \^A^''-^\\\Tkek\\} for fc > 1. (6.1) 

Then P|l > \\Tk\\ > A'-''\ Clearly, A'-'"'^ is monotonically increasing and is thus an 
improving estimate for as k increases. By the property of QLP decomposition in 
(4.2) and (5.2), we could easily extend (6.1) to include the largest diagonal of Lk'- 

^(0) ^ 0^ j^ik) ^ ^ax{^(^-i), lllfcefcll, 7^22, T^^, 17^11 A: > 1. (6.2) 

Some other schemes inspired by Larsen [31, section A. 6.1], Higham [25], and Chen 
and Demmel [8] follow. For the latter scheme, we use an implementation by Kaustuv 
[28] for estimating the norms of the rows of A. 

1. [31] ll Tfclli > llTfc ll 

2. [31] ^\\n^J\\\i>\\Tk\\ 

3. [31] ||Tj|| < llTfell for small j = 5 or 20 

4. [25] Matlab function NORMEST(A), which is based on the power method 

5. [8] max, ||/ij||/-ym, where hf is the ith row of AZ, each column of Z e IR"^™ 
is a random vector of ±l's, and m is a small integer (e.g., m — 10). 

Figure 6.1 plots estimates of \\A\\ for 12 matrices from the Florida sparse matrix 
collection [54] whose sizes n vary from 25 to 3002. In particular, scheme 3 above with 
j — 20 gives significantly more accurate estimates than other schemes for the 12 matri- 
ces we tried. However, the choice of j is not always clear and the scheme adds a little 
to the cost of MINRES-QLP. Hence we propose incorporating it into MINRES-QLP 
(or other Lanzcos-based iterative methods) if very accurate \\A\\ is needed. Otherwise 
(6.2) uses quantities readily available from MINRES-QLP and gives us satisfactory 
estimates for the order of ||A||. 

6.4. Matrix condition numbers. We again apply the property of the QLP 
decomposition in (4.2) and (5.2) to estimate ^(Ifc), which is a lower bound for k{A): 

7min ^ min{7i, 7^ 7min ^ min{7,„in, -fl^l^^ 7^-1, l7i^^ 1} for ^ > 3, 
k(") = 1, K('=) = niax|K('=-i),— 1 for/c>l. (6.3) 

I 7min J 

6.5. Solution norms. We derive a recurrence relation for \\xk\\ whose cost is as 
low as computing the norm of a 3- or 4- vector. 

Since \\xk\\ = ||VfePfcWfe|| = ||wfe||, we can estimate ||a;fe|| by computing Xk = \\uk\\- 
However, the last two elements of Uk change in Uk+i (and a new element iik+i is 
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□ |max.{l|r.€.|l,lT.I}- 
■■<>■■ ||-4|||/||.4|| 

■■■*■■■ yiK'^nll,- ||/1|||/||4|| 

-■x-'-l "dT,)- l|A|||/||,4|| 
-C^ -I MT20)- ||/1|||/||.4|| 

I K0HMEST(/1)- ||,4|||/P|| 

B Clieu Rud Demmel 



1 2 3 4 5 6 7 



9 10 11 12 



Fig. 6.1. Relative errors in different estimates of \\A\\. This figure can be reproduced by 
testminresQLPNormAS. 



added). We therefore maintain Xfc-2 by updating it and then using it according to 

xf22 = lllxi'ia t^fUl Xk = 4-1 cf. (5.8) and (5.9). 

(2) 

Thus Xk-2 increases monotonicaUy but we cannot guarantee that and its recurred 
estimate Xk are increasing, and indeed they are not in some examples (see Figure 8.1). 

6.6. Projection norms. Sometimes the projection of the right-hand side vector 
b onto ICk{A,b) is required (for example, see [46]). A simple recurrence relation 
is Luf. = ||Aa;j.||^ = ^t-i + '''k ^^'^ derive it in the same way as shown in 

Lemma 3.3. With loq = we have uJk = \\Axk\\ — \\[u!k-i Tk]\\. 

7. Preconditioned MINRES and MINRES-QLP. It is often asked: How can 
we construct a preconditioner for a linear system solver so that the same problem is 
solved with fewer iterations? Previous work on preconditioning the symmetric solvers 
CG, SYMMLQ, or MINRES includes [43, 37, 17, 12, 14, 35, 42, 34, 20, 2, 52]. 

We have the same question for singular symmetric systems Ax ~ b. Two-sided 
preconditioning is needed to preserve symmetry. We can still solve compatible sys- 
tems, but we will no longer obtain the minimum-length solution. For incompatible 
systems, preconditioning alters the "least squares" norm. To avoid this difficulty we 
must work with larger equivalent systems that are compatible. We consider each 
case in turn, using a positive-definite preconditioner M — CC^ with MINRES and 
MINRES-QLP to solve symmetric compatible systems Ax = b. Implicitly, we are solv- 
ing equivalent symmetric systems C^^AC~^y = C~^b, where C^x = y. As usual, it is 
possible to work with M itself, so without loss of generality we can assume C = Af 2 . 

7.1. Derivation. We derive preconditioned MINRES for compatible Ax — b 
by applying MINRES to the equivalent problem Ax = b, where A = M~2AM~2^ 
b = M~2b, and x = M~2x. 

7.1.1. Preconditioned Lanczos process. Let Vk denote the Lanczos vectors 
of /C(A, b). With vq — and f3iVi = 6, for fc = 1, 2, . . . we define 

Zk=l3kMhk, qk= PkM-^Vk, so that Mqk = Zk- (7.1) 
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Then 13k = ||^fcffe|| = \\M ^Zk\\ = WzkWm-^ = IkfelU-/ = yli^k, where the square root 
is well defined because M is positive definite, and the Lanczos iteration is 

Pk = Avk = M-'^AM-'^Vk = Af~5Agfe//3fc, 
"fe = vlpk = qlAqk/Pl, 
Pk+iVk+i = M^'^AM^'^Vk - akVk - PkVk-i- 

Multiplying the last equation by we get 

Zk+i = Pk+iM^Vk+i = AM^hk - akMKk - PkMKk-i 

1 ■ "fe /3fc 

The last expression involving consecutive Zj^s replaces the three-term recurrence in 
Dj 's. In addition, we need to solve a linear system Mqk = Zk (7.1) each iteration. 

7.1.2. Preconditioned MINRES. From (3.4) and (3.8) we have the following 
recurrence for the kth column of Dk = VkRk^ and Xk'- 

dk = {vk - 5f ■'dfc-i - ^kdk-2)hf\ Xk = Xk-i + Tkdk. 

Multiplying the above two equations by M^s on the left and defining dk = M~^dk, 
we can update the solution of our original problem by 



dk = (y^-^qk-5'i^''dk-^i~ekdk~2^/lf\ Xk = M Xk ^ Xk-i + Tkdk- 

We hst the algorithm in [10, Table 3.4]. 

7.1.3. Preconditioned MINRES-QLP. A preconditioned MINRES-QLP can 
be derived very similarly. The additional work is to apply right reflectors Pk to Rk, 
and the new subproblem bases are Wk = VkPk, with Xk — WkUk- Multiplying the 
new basis and solution estimate by M~2 on the left, we obtain 

Wk = M-'^Wk ^ M-'^VkPk, 

— -- — - (2) (2) -(3) -(2) 

Xk^M ^Xk^M ^T^feMfc = VKfeMfe = a;^^2 + A^fe-i^fc-i + A^feK • 

Algorithm 1 lists all steps. Note that Wk is written as Wk for all relevant k. Also, the 
output X solves Ax ~ b but the other outputs are associated with Ax ~ b. 

Remark. The requirement of positive-definite preconditioners M in MINRES and 
MINRES-QLP may seem unnatural for a problem with indefinite A because we cannot 
achieve M~2AM~2 « /. However, as shown in [17], we can achieve AI^^AM^^ « 
using an approximate block-LDL""" factorization A « LDLF to get M = 
L\D\L'^, where D is indefinite with blocks of order 1 and 2, and |D| has the same 
eigensystem as D except negative eigenvalues are changed in sign. 

SQMR [15] without preconditioning is analytically equivalent to MINRES. Unlike 
MINRES, SQMR can work directly with an indefinite preconditioner (such as block- 
LDL"""). However, in finite precision, SQMR needs "look-ahead" to prevent numerical 
breakdown. 
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Algorithm 1: Preconditioned MINRES-QLP to solve {A - al)x « b. 
input: A, b, a, M 

1 zo = 0, zi = b, Solve Mqi = zi, /3i = y/Wq\ [Initialize] 

2 Wo = W-\ = 0, X-2 = X-1 = xo — 

3 Co,l=Co,2=Co,3=-l, SO,l='SO,2=S0,3=0, </)0=A, = = X-2 = X- 1 = XO = 

4 Si = 7_i = 70 = n-i =Vo ='01= '^-i = '&0 = •&! = //-I = /^o = 0, ^ = 0, k = 1 

5 fc = 

6 while no stopping condition is satisfied do 
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fc A; + 1 

Pk = Aqk — aqk, ak = -kiquPk [Preconditioned Lanczos] 

"fc 

Zk+l = j-^Pk - f-^Zk - -^^Zk-l 

Solve Mqk+1 = Zk+1, Pk+i = "^ql+i^k+i 

if A; = 1 then pfc = II [afc ;9fc+i]|| else = ||[/3fc ak l3k+i]\\ 

(2) 

^k ''k-i.i^k + Sk-i.i<^k [Previous left reflection...] 

7fc = Sfe_i,i(5fe — Cfe_i,iQfc [on middle two entries of TkCk . . .] 

€fc+i = Sk-i,ifik+i [produces first two entries in Tfe+iefc+i] 

Sk+i = — Cfe_i,i/3fc+i 

Cfci, Sfci, 7^^' SymOrtiio(7fe, [Current left reflection] 
Cfc2,Sfe2,7^'!^2 SymOrtlio(7^i^2'«fe) [First right reflection] 

x(3) a j>(2) (3) (2) (2) 

Ok = Sk2l!fk-1 - Ck20l ' , 7fc — -Cfc27fc 5 r)k=Sk2Tk 

■&fll = Cfe2l?fe-1 + 3k2&^k'' 

Cfc3,Sfe3,7fcl'i ^ SymOrtiio(7^'i\,<5^^') [Second right reflection...] 

a (») (4) (3) j_ r(3)T 

Vk = Sfe37fe , 7fc = — Cfc37^ [to zero out 6). '] 

Tfc — Cki4>k-i [Last element of ffc] 

4>k = Ski<t>k-i, iik-i = 4>k-i\\['yk <5fc+i]|| [Update ||rfc||, ||^rfe_i||] 

if A; = 1 then 7inin = 7i else 7min •<- min {7min, 7fe-2' li-n l7fe'''|} 

= ma.c{^('=-i),p,,7il^2,7£\, 17^1} [Update ||^||] 

iOk = \\[(^k-i -Tfelll, K -s- ^^*'V7mm [Update IIAxfell, k(^)] 

Wfe = -(Cfe2//3fe)gfe + Sfe2wfi2 [Update Wk-2> Wk-i, Wkl 

wi-2 = iSk2/l3k)qk+Ck2W^k-2 

•CI ^ n 4.1. (2) (2) (3) (2) , 

if A; > 2 then w}, = SfcsWj..^ — CksWk, w^Ji = Cfc3™fe_i + SkaWk 
if A; > 2 then ^^^ij = {Tk-2 - 'nk-2fJ^^k-i ~ '^'=-2Mfe^^3)/7fc-2 [Update ^^-2] 
if A: > 1 then yU*.^2i = (jk-i - nk-iH^k-z ~ "^^k-l^J'^k-2) hk'-i [Update /x^-i] 
if 7fc^' 7^ then /xfc = (rjb - Jjfc/ii^ij ~ '^kljfkli)h^k^ else fik = [Compute /xfc] 
3^1-2 = a;^^23 + /if22w;r22 [Update Xk-2'] 

(2) , (2) (3) , (2) rr. 4- n 

a;*; = xl^2 + Mfc-i™fc-i + MfcW^ [Compute Xki 

& = \\[& ^kUW [Update ||a;fe_2||] 

Xfc = l|[Xfc-2 Mfc-i [Compute ||a;/c||] 



37a; = XA:, 4> = <(>k, V' = ^fc||[7fc+i '5fc+2]||, X = Xk, A = A'-'^\ oj -- 

output: X,(j),1p,X:A,K,U) 

38 

[c, s •<— SymOrtho(o, 6) is a stable form for computing r = \/a? + 6^, c ^ 
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7.2. Preconditioning singular Ax — b. For singular compatible systems, MIN- 
RES and MINRES-QLP find the minimum-length solution (see Theorems 3.1 and 5.1). 
If M is nonsingular, the preconditioned system is also compatible and the solvers 
return its minimum-length solution. The unpreconditioned solution solves Ax ~ b, 



Then rank(A) = 3 and Ax = b 



but is not necessarily a minimum-len 


rth solution. 




■ 1 1 0' 




'6' 


Example 7.1. Let A = 


1110 
10 1 
10 


and b — 


9 
6 
3 



is a singular compatible system. The minimum-length solution is x 



t - 



2 4 3 2 



By 



binormalization [33] we construct the matrix D = diag([o. 84201 o. 81228 0.30957 3.2303]). 
The minimum-length solution of the diagonally preconditioned problem DADy — Db 

is y^ = [3.5739 3.6819 9.6909 0.93156 J"^' Then X = Dy^ ~ [3.0092 2.9908 3.0000 3.0092]"^ is a 

solution of Ax = b, but x ^ x^ . 

7.3. Preconditioning singular Ax b. We propose the following techniques 
for obtaining minimum-residual solutions of singular incompatible problems. In each 
case we use an equivalent but larger compatible system to which MINRES may be 
applied. Even if the larger system is singular, Theorem 3.1 shows that the minimum- 
length solution of the larger system will be obtained. The required x will be part of 
this solution. Preconditioning still gives a minimum-residual solution of Ax ~ b, and 
in some cases x will be x^. If the systems are ill-conditioned, it will be safer and 
more efficient to apply MINRES-QLP to the original incompatible system. However, 
preconditioning will give an x that is "minimum length" in a different norm. 

7.3.1. Augmented system. When A is singular, so is the augmented system 



(7.2) 



but it is always compatible. Preconditioning with symmetric positive-definite M gives 
us a solution [^] in which r is unique, but x may not be x^ . 

7.3.2. A giant KKT system. Problem (1.1) is equivalent to min^, x x^x subject 
to (7.2), which is an equality-contrained convex quadratic program. The correspond- 
ing KKT system [36, section 16.1] is both symmetric and compatible: 
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.0. 



(7.3) 



Although this is still a singular system, the upper-left 3x3 block-submatrix is nonsin- 
gular and therefore r, x, and y are unique and a preconditioner applied to the KKT 
system would give x as the minimum-length solution of our original problem. 

7.3.3. Regularization. If the rank of a given matrix A is ill-determined, we 

may want to solve the regularized problem [13, 22] with parameter 5 > 0: 



mm 



A 
51 



(7.4) 



The matrix [^] has full rank and is always better conditioned than A. LSQR [40, 41] 
may be applied, and its iterates Xk will reduce + ^^ll^Cfcll^ monotonically. Alter- 

natively, we could transform (7.4) into the following symmetric compatible systems 
and apply MINRES or MINRES-QLP. They tend to reduce \\Ark-5'^Xk\\ monotonically. 
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{A^ + 6^I)x = Ab. 



Normal equation: 
Augmented system: 

A two-layered problem: If we eliminate v from the system 



(7.5) 
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(7.6) 



we obtain (7.5). Thus x is also a solution of our regularized problem (7.4). 
This is equivalent to the two-layered formulation (4.3) in Bobrovnikova and 
Vavasis [5] (with Ai = A, A2 ^ Di ^ D2 = I, &i = 6, 62 = 0, 5i = 1, 
= (5^). A key property is that a; — as J — ?> 0. 
A KKT-like system: If we define y — —Av and r — b — Ax — 5^y, then we can show 
(by eliminating r and y from the following system) that x in 
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(7.7) 



is also a solution of (7.6) and thus of (7.4). The upper- left 3x3 block- 
submatrix of (7.7) is nonsingular, and the correct limiting behavior occurs: 
a; — > as (5 — > 0. In fact, (7.7) reduces to (7.3). 

7.4. General preconditioners. The construction of preconditioners is usually 
problem-dependent. If not much is known about the structure of A, we can only 
consider general methods such as diagonal preconditioning and incomplete Cholesky 
factorization. These methods require access to the nonzero elements of A. (They are 
not applicable if A exists only as an operator for returning the product Ax.) 

For a comprehensive survey of preconditioning techniques, see Benzi [3]. We 
discuss a few methods for symmetric A that also require access to the nonzero Aij. 

7.4.1. Diagonal preconditioning. If A has entries that are very different in 
magnitude, diagonal scaling might improve its condition. When A is diagonally dom- 
inant and nonsingular, we can define D = diag((ii, . . . , (i„) with dj — \/\Ajj\^/'^ . 
Instead of solving Ax = 5, we solve DADy = Db, where DAD is still diagonally 
dominant and nonsingular with all entries < 1 in magnitude, and x — Dy. 

More generally, if A is not diagonally dominant and possibly singular, we can 
safeguard division-by-zero errors by choosing a parameter S > and defining 



dj{S) — 1/ max{(5, \/ \ Ajj \ , max|ylij|}, j = 1, . . . ,n. 



(7.; 



Example 7.2. 
1. IfA - 



1 10" 

-s 



10"° 1 10" 
10* 



m (7.8). Then DAD 



then k{A) w 10"*. Let 5 ^ I, D ^ 
and k{DAD) « 1. 



10" 



10" 



-1 10" 
-1-10 in-^l 
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2. A 



10" 
10" 



10" 
10" 
10" 



contains mostly very small entries, and k{A) 10 



10 



Let S = 10-^ and D 



10^ 



10^ 



10° 



10" 



Then DAD 



10" 



10" 

1 

10^ 



10^ 




and k(DAD) w 10^. (The choice of S makes a critical difference in this case: 
with 5 — 1, we have D — I .) 

7.4.2. Binormalization (BIN). Livne and Golub [33] scale a symmetric ma- 
trix by a series of k diagonal matrices on both sides until all rows and columns of the 
scaled matrix have unit 2-norm: DAD = D^ ■ ■ ■ DiADi ■ ■ ■ Dk- See also Bradley [6]. 



Example 7.3. If A ^ 
sweep of BIN, we obtain D ~ 



10" 



10 



then K,{A) 
diag(8.1e-3,6.6e-5,1.5), DAD 



12 



With just one 



1 

10"** 10 
10"* 

" 6.5e-l 5.3e-l 
5.3e-l 1 
1 

and k(DAD) « 2.6 even though the rows and columns have not converged to one 
in the two-norm. In contrast, diagonal scaling (7.8) defined by 5 = 1 and D — 



diag(l, 10" 



) reduces the condition number to approximately 10 . 



7.4.3. Incomplete Cholesky factorization. For a sparse symmetric positive 
definite matrix A, we could compute a preconditioner by the incomplete Cholesky 
factorization that preserves the sparsity pattern of A. This is known as ICO in the 
literature. Often there exists a permutation P such that the ICO factor of PAP^ is 
more sparse than that of A. 

When A is semidefinite or indefinite, ICO may not exist, but a simple variant that 
may work is the incomplete Cholesky-infinity factorization [55, section 5]. 

8. Numerical experiments. We compare the computed results of MINRES- 
QLP and various other Krylov subspace methods to solutions computed directly by 
the eigenvalue decomposition (EVD) and the truncated eigenvalue decompositions 
(TEVD) of A. For TEVD we have 



1 



|A,|>t||A||e 



ll^ll = max|A^|, KtiA) = 



max I Xi I 
min I Xi 

|A,|>t||A||e 



with parameter t > 0. Often t is set to 1, and sometimes to a moderate number 
such as 10 or 100; it defines a cut-off point relative to the largest eigenvalue of A. 
For example, if most eigenvalues are of order 1 in magnitude and the rest are of 
order [j^lje « 10^^^, we expect TEVD to work better when the small eigenvalues are 
excluded, while EVD (with t = 0) could return an exploding solution. 

In the tables of results, Matlab MINRES and Matlab SYMMLQ are Matlab's 
implementation of MINRES and SYMMLQ respectively. They incorporate local re- 
orthogonalization of the Lanczos vector V2, which could enhance the accuracy of the 
computations if b is close to an eigenvector of A [32] : 

Second Lanczos iteration: Pivi = 6, and (72 = P2V2 — Avi — aiVi 
Local reorthogonalization: 172 ^ 92 ^ {viq2)vi. 

Lacking the correct stopping condition for singular problems, Matlab SYMMLQ 
works more than necessary and then selects the smallest residual norm from all com- 
puted iterates; it would sometimes report that the method did not converge although 
the selected estimate appeared to be reasonably accurate. 
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MINRES SOL and SYMMLQ SOL are implementations based on [39]. MINRES+ 
and SYMMLQ+ are the same but with additional stopping conditions for singular 
incompatible systems (see Lemma 3.3 and [10, Proposition 2.12]). 

The computations in this section were performed on a Windows XP machine with 
a 3.2GHz Intel Pentium D Processor 940 and 3GB RAM (e « lO^^*^) . Tests were 
performed with each solver on five types of problem: 

1. symmetric, nonsingular linear systems 

2. symmetric, singular linear systems 

3. mildly incompatible symmetric (and singular) systems (meaning ||r|| is rather 
small with respect to ||6||) 

4. symmetric (and singular) LS problems 

5. Hermitian systems. 

We present a few examples that illustrate the key features of MINRES-QLP. For 
a larger set of tests and results, such as applying MINRES-QLP and other Krylov 
methods to Hermitian systems with preconditioning, we refer to [10, Chapter 4]. 

For a compatible system, we generate a random vector b that is in the range 
of the test matrix (6 = Ay, yi ~ i.i.d. J7(0, 1), i.e., yi,...,j/„ are independent and 
identically distributed random variables, whose values are drawn from the standard 
uniform distribution with support [0, 1]). For an LS problem, we generate a random 
h that is not in the range of the test matrix (bi ^ i.i.d. L/(0, 1) often suffices). 

If A is Hermitian, then v^Av is real for all complex vectors v. Numerically (in 
double precision), = vjf Avk is likely to have small imaginary parts in the first few 
Lanczos iterations and snowball to have large imaginary parts in later iterations. This 
would result in a poor estimation of ||Tfe||i? or H^Hf, and unnecessary errors in the 
Lanczos vectors. Thus we made sure to typecast ak = rea.\{vlf Avk) in MINRES-QLP 
and MINRES SOL. 

We could say from the results that the Lanczos-based methods have built-in 
regularization features [29], often matching the TEVD solutions very well. 

8.1. A Laplacian system Ax « b (almost compatible). Our first example 
involves a singular indefinite Laplacian matrix A of order n — 400. It is block- 
tridiagonal with each block being a tridiagonal matrix T of order iV = 20 with all 
nonzeros equal to 1: 



A = 



T T 
T T 



■■ T 

T T 



T = 



1 1 
1 1 



(8.1) 



NxN 



Matlab's eig(A) reports the following data: 205 positive eigenvalues in the inter- 
val [6.1e— 2, 8.87], 39 almost-zero eigenvalues in [— 2.18e— 15, 3.71e— 15], 156 negative 
eigenvalues in [—2.91, — 6.65e— 2], numerical rank — 361. 

We used a right-hand side with a small incompatible component: b = Ay+ 10~^z 
with yi and Zi ~ i.i.d. U(0, 1). Results arc summarized in Table 8.1. In the column 
labeled "C?", the value "Y" denotes that the associated algorithm in the row has 
converged to the desired NRBE tolerances within maxit iterations (cf. Table 6.1); 
otherwise, we have values "N" and "N?", where "N?" indicates that the algorithm 
could have converged if more relaxed stopping conditions were used. The column 
"Aw" shows the total number of matrix-vector products, and column "a;(l)" lists the 
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Table 8.1 

Finite element problem Ax b with b almost compatible. Laplacian on a 20 X 20 grid, n = 400, 
maxit = 1200, shift = 0, tol = l.Oe— 15, maxnorm = 100, maxcond = lel5, ||6|| = 87. To reproduce 
this example, run test_minresqlp_eg7_l (2^.) ■ 
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first element of the final solution estimate x for each algorithm. For GMRES, the 
integer in parentheses is the value of the restart parameter. 

MINRES SOL gives a larger solution than MINRES-QLP. This example has a 
residual norm of about 1.7 x 10~^, so it is not clear whether to classify it as a linear 
system or an LS problem. To the credit of Matlab SYMMLQ, it thinks the system 
is linear and returns a good solution. For MINRES-QLP, the first 410 iterations are in 
standard "MINRES mode" , with a transfer to "MINRES-QLP mode" for the last 202 
iterations. LSQR [40, 41] converges to the minimum-length solution but with more 
than twice the number of iterations of MINRES-QLP. The other solvers fall short in 
some way. 

8.2. A Laplacian LS problem min||Aa; — 6||. This example uses the same 
Laplacian matrix A (8.1) but with a clearly incompatible 6 = 10 x rand(n, 1), i.e., 
bi ^ i.i.d. U{0,10). The residual norm is about 17. Results are summarized in 
Table 8.2. MINRES gives an LS solution, while MINRES-QLP is the only solver that 
matches the solution of TEVD. The other solvers do not perform satisfactorily. 

8.3. Regularizing efTect of MINRES-QLP. This example illustrates the reg- 
ularizing effect of MINRES-QLP with the stopping condition Xk < maxxnorm. For 
A: > 18 in Figure 8.1, we observe the following values: 

3.87e-ll 1.38 X 10^] II 1.38 X 10^ 

8.00e-10 -1.52 x 10^] || 1.52 x 10^ 

1.62e-10 -1.62 x 10^] || = 1.62 x 10® > maxxnorm = lO". 

Because the last value exceeds maxxnorm, MINRES-QLP regards the last diagonal 
element of Lk as a singular value to be ignored (in the spirit of truncated SVD 
solutions). It discards the last element of M20 and updates 

X2o^||[2.51 1.62e-10 0]||=2.51. 



Xi8 = II [2.51 
Xi9-||[2.51 - 
X20 = II [2.51 
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Table 8.2 

Finite element problem min||Ax — Laplacian on a 20 X 20 grid, n = 400, maxit = 500, 
shift = 0, tol = l.Oe— 14, maxnorm = le4, maxcond = lel4, = 120. To reproduce this example, 
run test_minresqlp_eg7_l (25) . 
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The full truncation strategy used in the implementation is justified by the fact 
that Xk = WkUk with Wk orthogonal. When becomes large, the last element of 
Uk is treated as zero. If ||a;/c|| is still large, the second-to-last element of Uk is treated 
as zero. If HxfcH is still large, the third-to-last element of Uk is treated as zero. 

8.4. Effects of rounding errors in MINRES-QLP. The recurred residual 
norms cf)^^ in MINRES usually approximate the directly computed ones very 
well until ||r^^|| becomes small. We observe that 4>^' continues to decrease in the last 
few iterations, even though has become stagnant. This is desirable in the sense 

that the stopping rule will cause termination, although the final solution is not as 
accurate as predicted. 

We present similar plots of MINRES-QLP in the following examples, with the cor- 
responding quantities as 0^ and ||?'^||. We observe that except in very ill-conditioned 
LS problems, 0^ approximates very closely. 

Figure 8.2 illustrates four singular compatible linear systems. 

Figure 8.3 illustrates four singular LS problems. 

9. Conclusion. MINRES constructs its fcth solution estimate from the recursion 
Xk — Dktk — Xk-i + Tkdk (3.4), where n separate triangular systems RJD'^ = Vj[ 
are solved to obtain the n elements of each direction di, . . . ,dk. (Only dk is obtained 
during iteration k, but it has n elements.) 

In contrast, MINRES-QLP constructs its /cth solution estimate using orthogonal 
steps: x2 = {VkPk)uk; see (5.3)-(5.9). Only one triangular system LkUk — QkiPiCi) 
is involved for each k. 

Thus MINRES-QLP overcomes the potential instability predicted by the MINRES 
authors [39] and analyzed by Sleijpen et al. [47]. The additional work and storage 
are moderate, and maximum efficiency is retained by transferring from MINRES to 
the MINRES-QLP iterates only when the estimated condition number of A exceeds a 
specified value. 
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Fig. 8.1. Recurred if>k « ||rfc||, i/ij, ki \\Ark\\, and \\xk\\ for MINRES and MINRES-QLP. The 
matrix A (ID 1177 from [54]) is positive semidefinite, n = 25, and b is random with \\b\\ ~ 1.7. 
Both solvers could have achieved essentially the TEVD solution of Ax ~b at iteration 11. However, 
the stringent tol = 10~^* on the recurred normwise relative backward errors (NRBE in Table 6.1) 
prevents them from stopping "in time" . MINRES ends with an exploding solution, yet MINRES- 
QLP brings it back to the TEVD solution at iteration 20. Left: <p^' and (recurred ||rj,|| of 
MINRES and MINRES-QLP) and their NRBE. Middle: ip^' and tp^ (recurred \\Ark\\) and their 
NRBE. Right: ||a;^'^|| (norms of solution estimates from MINRES) and (recurred \\xf;\\ from 
MINRES-QLP) with maxxnorm = 10*. This figure can be reproduced by test-minresqlp-fig7-l (2) . 



MINRES and MINRES-QLP are readily applicable to Hermitian matrices, once ak 
is typecast as a real scalar in finite-precision arithmetic. For both algorithms, we 
derived recurrence relations for ||Arfc|| and ||Aa;j,|| and used them to formulate new 
stopping conditions for singular problems. 

TEVD or TSVD are commonly known to use rank-fc approximations to A to find 
approximate solutions to min \ \Ax — b\\ that serve as a form of regularization. Krylov 
subspace methods also have regularization properties [23, 21, 29]. Since MINRES- 
QLP monitors more carefully the rank of T^,, which could be k or k — 1, we may say 
that regularization is a stronger feature in MINRES-QLP, as we have shown in our 
numerical examples. 

It is important to develop robust techniques for estimating an a priori bound for 
the solution norm since the MINRES-QLP approximations are not monotonic as is the 
case in CG and LSQR. Ideally, we would also like to determine a practical threshold 
associated with the stopping condition 7^^'' = in order to handle cases when 7^^'' is 
numerically small but not exactly zero. These are topics for future research. 

10. Software and reproducible research. Matlab 7.6, Fortran 77, and For- 
tran 90 implementations of MINRES with new stopping conditions ||Arfc|| < A|| ||rfc|| 
and \\Axk\\ < tol\\A\\\\xk\\, and a Matlab 7.6 implementation of MINRES-QLP are 
available from SOL [48] . 

Following the philosophy of reproducible computational research as advocated in 
[11, 9], for each figure and example in this paper we mention either the source or the 
specific Matlab command. Our Matlab scripts are available at SOL [48]. 
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Four symmetric positive semidefinite ill-conditioned systems. 



K i« 10', \\A\\ = 3, \\x\\ a in«, ||!,|| as lO'.M = Kl"" 
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Fig. 8,2. Solving Ax = b with semidefinite A similar to an example of Sleijpen et al. [47]. 
A = Q diag([05, ?7, 2r;, 2 : : 3])Q of dimension n = 797, nullity 5, and norm \\A\\ = 3, where 
Q = I — (2/n)ww'^ is a Householder matrix generated by v = [O5, 1, . . . , 1]^, w = v/\\v\\. These plots 
illustrate and compare the effect of rounding errors in MINRES and MINRES-QLP. 

The upper part of each plot shows the computed and recurred residual norms, and the lower part 
shows the computed and recurred normwise relative backward errors (NRBE, defined in Table 6.1). 
MINRES and MINRES-QLP terminate when the recurred NRBE is less than the given tol = lO"^''. 

Upper left: r) = 10~^ and thus k(A) 10*. Also b = e and therefore \\x\\ 2> \\b\\. The graphs 
of MINRES's directly computed residual norms \\r^^\\ and recurrently computed residual norms 
(f>jf start to differ at the level of 10~^ starting at iteration 21, while the values </<^ ||r^|| from 
MINRES-QLP decrease monotonically and stop near 10~® at iteration 26. 

Upper right: Again rj = 10~* but b = Ae. Thus \\x\\ = ||e|| = 0{||b|i). The MINRES graphs 
of ||r^'^|| and <j)^' start to differ when they reach a much smaller level of lO"-'" at iteration 30. The 
MINRES-QLP 4>^ 's are excellent approximations of (j>^ , with both reaching lO"-*^"^ at iteration 33. 

Lower left: r] = 10^^*^ and thus A is even more ill-conditioned than the matrix in the upper 
plots. Here b = e and \\x\\ is again exploding. MINRES ends with 10'^, which means no 

convergence, while MINRES-QLP reaches a residual norm o/10~*. 

Lower right: r) = 10"-"^ and b = Ae. The final MINRES residual norm ||r^^|| Ri 10~*, which 
is satisfactory but not as accurate as <p^^ claims at lO"-"-"^. MINRES-QLP again has (j)^ Ik^ll ~ 
10~^^ at iteration 37. 

This figure can be reproduced by DPtestSingJ.m. 
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Four (singular) symmetric least squares problems. 

: 10M|-4|1 =3,||.t|1 ~ 10", ||6|| ~ W\tol = 10-1 H~ 10', II All =3,||ii| =s 10", ||6|| 10', fo! = lO"" 





5 10 15 20 25 



10". \\A\\ =3, ||xi| » 10", Hill » 10', fo/ = 10-" 



15 

k 

Fig. 8.3. Solving Ax 
Qdiag([05,)7,2»?,2 




1 

789 



30 



= b with semidefinite A similar to an example of Sleijpen et al. [47]. 
:3])Q of dimension n = 797 with \\A\\ = 3, where Q = I — {2/n)ee'^ is a 



Householder matrix generated by e = [1, . . . , 1] . (We are not plotting the NRBE quantities because 



|rf;|| Si 6 throughout the iterations in this example.) 
Upper left: r) = 10~^ and thus cond(A) lO'^. Also b = e and therefore \\x\\ 3> The 
graphs of MINRES's directly computed ||Ar^|| and recurrently computed tp^^ , and also tp'^ |jAr^|| 
from MINRES-QLP, match very well throughout the iterations. 

Upper right: Here, rj = lO"** and A is more ill-conditioned than the last example (upper left). 
The final MINRES residual norm 11^?-^^ II is slightly larger than the final MINRES-QLP 

residual norm ^? Si \\Ar^\\. The MINRES-QLP V? are excellent approximations of ||Ar?||. 



' k 

10 



and cond(A) ^ 10''. MINRES's Vf and 



ArfW 



M'-f II 



fer starting at 



3, which means no convergence. MINRES-QLP reaches a 



Lower left: rj 

iteration 21. Eventually, \\ 
residual norm of = ||ArJ^|| = lO"'^. 

Lower right: rj = 10~®. MINRES performs even worse than in the last example (lower left). 
MINRES-QLP reaches a minimum ||Ar?|| Oi 10"^ but tol = 10~^ does not shut it down soon enough. 



The final ipj^ 



WAri 



10- 



The values of tp'j^ and ||Arj?|| differ only at iterations 27-2 



This figure can be reproduced by DPtestLSSingS . m. 
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Appendix A. Proof that Tg is nonsingular iff 6 e rangc(74) (section 2.1). 

If Ti is nonsingular, we have AViT^^ei = Vgei = /3f Conversely, if 6 € 
range(A), then range(Vf) C range(yl) and null(A) n range(y£) — {0}. We also know 
that Tank{Vi) = £ and rank(Ti;) = Tauk{VeTi) = Ta,uk{AVe) = rank(V£)-dim[null(A)n 
range(T/f)]; see [4, Fact 2.10.4 ii]. Thus rank(T^) — £ and so T( is nonsingular.) 
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